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Abstract. In the context of forecasting temperature and pressure fields in high-intensity 
focussed ultrasound, the accuracy of predictive models is critical for the safety and efficacy 
of treatment. In such fields inertial cavitation is often observed. Classically, estimations of 
cavitation thresholds have been based on the assumption that the incident wave at the surface 
of a bubble was the same as in the far-field, neglecting the effect of nonlinear wave propagation. 
By modelling the incident wave as a solution to Burgers' equation using weak shock theory, the 
effects of nonlinear wave propagation on inertial cavitation are investigated using both numerical 
and analytical techniques. From radius-time curves for a single bubble, it is observed that there is 
a reduction in the maximum size of a bubble undergoing inertial cavitation and that the inertial 
collapse occurs earlier in contrast with the classical case. Corresponding stability thresholds for a 
bubble whose initial radius is slightly below the critical Blake radius are calculated. Bifurcation 
diagrams and frequency-response curves are presented associated with the loss of stability. The 
consequences and physical implications of the results are discussed with respect to the classical 
results. 



1. Introduction 

For materials with a nonlinear stress-strain relationship, such as tissue [1], the point of maximum 
compression for a wave may propagate faster than the point of maximum rarefaction, leading 
to a distortion in the wave profile and the redistribution of energy from the fundamental 
harmonic frequency to higher harmonics. The pressures associated with therapeutic high 
intensity focused ultrasound may be high enough for the effects of nonlinear wave propagation to 
be significant [2] . As high frequency components are absorbed more easily than those with lower 
frequencies, nonlinear wave propagation contributes to increased absorption, enhanced heating 
and a subsequent shift in the focal point of targeted ultrasonic beam, potentially damaging 
healthy tissue. As well as providing greater predictive accuracy, knowledge of nonlinear wave 
propagation will enable increased accuracy in calibration using the received signal generated 
by bubble oscillations [3] and thus greater accuracy in treatment procedures. In this paper the 
implications of nonlinear wave propagation on inertial cavitation are investigated and stability 
criteria re-derived in an attempt to reconcile theory and experiment. 

The Rayleigh-Plesset equation [4] is a nonlinear equation which determines the size of a 
spherical bubble subject to a varying pressure field. Wave propagation from the far-field to 
the surface of the bubble is generally assumed to be linear, yet this does not correspond with 
experimental observations in the context of high intensity focused ultrasound (HIFU). Moss [5] 



attempted to incorporate the effects of boundary conditions and global compressibility /local 
incompressibility into the Rayleigh-Plesset equation but only considered linear wave propagation. 
The equation derived was identical to the classical Rayleigh-Plesset equation if the far-field 
pressure is replaced by an attenuated pressure at the bubble surface. In this paper the effect of 
distortion rather than attenuation of the wave profile and its effect on the stability of oscillations 
will be investigated. 

Lauterborn and co-workers [4, 6] and Smerka [7] showed experimentally and numerically that 
bubble oscillations undergo a sequence of period-doubling bifurcations leading to unstable quasi- 
periodic (chaotic) oscillations. The period doubling route to chaos occurs through a succession 
of saddle-node bifurcations of subharmonic periodic orbits. Homoclinic bifurcations are the limit 
of a countable sequence of subharmonic saddle-node bifurcations and thus provide an insight 
into the parameters at which unbounded growth may begin to occur. 

Mel'nikov's method provides a measure of the distance between the stable and unstable 
manifolds of a periodically perturbed system. If the manifolds intersect transversely once, 
they will intersect transversely infinitely many more times. The transverse intersections can 
be represented by Smale horseshoes, which through the Smale-Birkhoff theorem give an elegant 
description of sensitivity to initial conditions and the resulting chaotic oscillations [8] . The first 
application of Mel'nikov's method to cavitation was performed by Chang and Chen [9] on the 
effect of viscosity on the Hamiltonian structure. Harkin [10] also performs Mel'nikov analysis 
on bubbles whose initial radius is slightly smaller than Blake critical radius. Using matched 
perturbation analysis, Harkin derives a second order normal form for the Rayleigh-Plesset 
equation. The normal form is a damped driven oscillator. An escape velocity, like the static Blake 
criterion, provides an upper bound for when unbounded growth will occur, whereas Mel'nikov's 
method provides a lower bound for when the transition to chaos and unbounded growth may 
occur. The fate of bubbles whose initial conditions lie in the intermediate region between 
the Mel'nikov and Blake thresholds can be computed by a transport-type processes [11]. A 
Bernoulli shift map on two symbols has already been constructed numerically from a bifurcation 
diagram for the full Rayleigh-Plesset equation [12] without explicit inference to sensitivity to 
initial conditions. In each application of Mel'nikov method the incident pressure wave was 
sinusoidal [7,9-13]. In this paper the analysis is performed for nonlinear waves. 

The outline of this paper is as follows, in section [2] nonlinear wave propagation is considered 
and wave profiles derived. Then in section [3] the Rayleigh-Plesset equation is introduced and the 
effects of nonlinear wave propagation investigated. In section 0] Mel'nikov analysis is performed 
for nonlinear wave propagation, providing an improved measure of the values at which quasi- 
periodic oscillation and unbounded growth may occur. Finally, conclusions and implications are 
discussed in section 

2. Nonlinear Wave Propagation 

There is no universally accepted system of partial differential equations for the modelling 
of ultrasound propagation in biological tissue. Perhaps the best known is the Khokhlov- 
Zabolotskaya-Kuznetsov (KZK) equation. The KZK equation is a parabolic wave equation which 
includes the effects of diffraction, absorption and nonlinearity of the directed beams [3]. The KZK 
equation for an axi-symmetric beam which propagates in the r direction is written in terms of 
the acoustic pressure p (r, t) as 
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where cq is the initial wave speed, p is the density of the medium, (3 is the standard nonlinearity 
coefficient given by (3 = 1 + B /2A, where B/A is the standard nonlinearity parameter, D is the 



sound diffusivity and t' = t — r/c is the retarded time variable. The Laplacian is taken with 
respect to the transverse coordinates. The sound diffusivity is given by 



D = - (fi b + 4/i s /3 + k (1/c - l/cp)) 

where /i& and /x s are the bulk and shear viscosity respectively, 1 < k < 7 the polytropic exponent, 
and c v and c p are the specific heats at constant volume and pressure respectively. If re = 1 the 
system is isothermal, if k = 7 = c p /c v the system is adiabatic. In the study of cavitation it will 
be assumed that the distance from the surface of the bubble to the shock front is a fixed length. 
The equation ([1]) can be integrated to give 

Discarding the effects of diffraction Burgers', or alternatively the Burgers-Hopf, equation is 
recovered 
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If D = the fluid is inviscid, and the so-called lossless Burgers equation is recovered. 

From a given driving pressure, f (t), the incident pressure wave, pit), can be expressed 
using weak shock theory. The location of a shock is determined by the Rankine-Hugoniot 
relation defining the conservation of flux. The areas enclosed by the multi-valued solution to 
the left and right of the shock are equal. Thus, by this symmetry, the shock is positioned at the 
zero of the linear solution. For a sinusoidal driving pressure of magnitude P and frequency u, 
i.e. / (t) = Psin (cot), a Fourier expansion of the solution to the lossless Burgers equation yields 
the Bessel-Fubini solution 
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where to n = nu>, J n is an n th -order Bessel functions of the first kind and r s is the normalised 
distance to a shock given by 

r pel 
r s = — where r c = —— 
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is the location of the shock. 

Beyond the shock, weak shock theory can once again be employed to find a solution, however, 
the resulting analytical solution takes an integral form. An asymptotic solution, valid for r s > 3 
is 



n=l 

which can be expressed in the time domain as 
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for - 7r < ut' < 0, 
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for < ut' < 7r. 



(6) 



A general solution to the Burgers'-Hopf equation for a sinusoidal driving pressure, derived 
by transforming the nonlinear equation into a linear diffusion equation via the Hopf-Cole 
transformation, is given by 
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where J n is an n th -order modified Bessel functions of the first kind and V is the Gol'dberg number 
which relates the importance of nonlinear effects against dissipative effects by 

r= 2/3P 
Dpuj 

If r ^> 1 then the relative effects of nonlinearity are strong, if T < 1 the relative effects of 
nonlinearity are weak. 

Unlike for strong nonlinearity, i.e. F> 1, the solution converges slowly but far from the 
shock, i.e. r s > 3, asymptotic analysis yields the more usable formulation 

_ 2P ^ sin (uj n t') 
P ~ r ^sinh(n(l + r s )/r) 1 j 

called the Fay solution. The Fay solution can be expressed as the periodic function 

In the lossless limit as viscosity vanishes, that is as T — > oo, the Fay solution recovers the Fourier 
expression for a sawtooth function ([6j). The Fubini (j3|) and Fay ([8]) solutions may at first seem 
incompatible but each holds in a different region of the flow; the Fubini solution close to the 
source as shocks begin to form and the Fay solution far from the source as shocks begin to 
decay. Blackstock [14] shows that the true solution to the lossless Burgers equation is simply the 
sum of the two solutions (jl]) and ([6]). In the near- field the Fubini solution is dominant, in the 
far-field the sawtooth solution is dominant. In [15] exact solutions of Burgers' equation from the 
Cole-Hopf transformation were computed numerically and contrasted against the solutions and 
showing good agreement between the sets of solutions unless in the immediate neighbourhood 
of the shock. Furthermore, the formation of a shock will result in a negative pressure and hence 
is a prime site for the nucleation of cavities. Thus the forthcoming analysis will be performed 
either before a shock or significantly after, so that only the effect of nonlinear wave propagation 
on pre-existing cavities is studied. 

Figure [T] contrasts the nonlinear wave profile and the linear wave profile illustrating the 
distortion due to nonlinear propagation of waves with equal amplitude. Throughout this paper 
comparisons between linear and nonlinear waves of equal amplitude are studied. 

3. Rayleigh-Plesset Equation 

The Rayleigh-Plesset equation is an ordinary differential equation which models the oscillations 
of a spherical bubble of radius R 

rr+^rA =PgiR)+Pv - PoD+pi t) + '^ + ^£ (10) 




Figure 1. Profiles of incident waves for nonlinear wave propagation. The nonlinear wave is 
modeled by the first twenty terms of the Bessel-Fubini solution. 



where p g is the internal pressure of the gas in the bubble given by the hardcore van der Waals 
relationship 

»W0) = (».-*i|)(|^)" (ID 

with Rq the equilibrium radius, a the van der Waals hard-core radius, p v the vapour pressure, p^ 
is the ambient pressure, a is the surface tension and \i viscosity. The gas is assumed to be ideal 
as the internal pressure is a function of the bubble radius only. The forcing pressure will take 
the form 

oo 

P (t) = ^ P n sin (u n (t + t )) 

71=1 

where P n are the Fourier terms of a solution to a nonlinear wave equation, such as those given 
by the Bessel-Fubini @ or the Fay ([H]) solutions. The Rayleigh-Plesset equation can be derived 
by balancing the energy supplied to the bubble by the incident pressure and the surround fluid 
and the kinetic energy of the bubble oscillations [16]. 

Figure contrasts the effects of nonlinear and linear wave propagation. It is clear that for 
nonlinear wave propagation inertial cavitation occurs before and at a smaller maximum radius 
than for linearly propagated waves. This has two significant effects, firstly as the collapse occurs 
at a small maximum radius there is a diminished chance of shape instability [17] and secondly, 
as the collapse occurs earlier more after-bounces can occurs so that the bubble returns to a 
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Figure 2. Contrasting profiles of inertial cavitation for linear (red) and nonlinear (blue) wave 
propagation modeled by the first ten terms of the Bessel-Fubini solution. Rq = 0.6/im, subject to 
a wave with frequency uj = 2ir ■ 26kHz, amplitude of driving pressure P = 1.36atm, normalised 
distance to the shock r s = 1/20 ambient pressure p^ = latm, surface tension a = 0.073kgm~ 2 , 
viscosity fi = 10 _3 kgm s _1 , density p = lOOOkgm" 3 , polytropic exponent k = 5/3 and hardcore 
radius a = Rq/8.85. 



stable initial radius before the next cycle. Note that the period of the after-bounces is the same 
for the nonlinear and linear wave since the after-bounces occur at the natural frequency of the 
Rayleigh-Plesset equation which is independent of the applied pressure. 

4. Mel'nikov Analysis 

Considering an unforced bubble and, discarding the effects of viscosity and the derivative of the 
gas pressure, let R (t) = Rq (1 + x). To first order, the nondimensional Rayleigh-Plesset equation 
yields the simple harmonic equation 

. + ^-0 where * -1 ( 3 " ~_ $ *> + 2. - . (12) 

The natural frequency oj is used to nondimensionalise time by r = u t so that the new time-like 
variable will be a function of the perturbation parameter. Note that when a = and k = 1 then 
the frequency given by Harkin [10, Eq. (12)] is recovered. In this section Harkin's analysis is 
followed as it gives an analytical expression for the transverse intersections which allows for 
comparison between linear and nonlinear wave propagation. 

The static Blake threshold pressure is the point at which the internal pressure p v + p g is 
equal to the external pressure p\ + 2a /R, thus for internal pressures larger than this threshold 
unbounded growth will occur. The equilibrium pressure exerted on the bubble surface by the 
liquid, pi, is given by 

2a . , 

Pl=Pg+Pv--j^- (13) 

Now perturb the equilibrium radius by R (t) = Rq (1 + (t)) with e a small parameter given 



by 

e = 2 (1 - R /R c ) 

where R c is the Blake critical radius, found as the stationary solutions to the unforced Rayleigh- 
Plesset equation 

2a where G = [poo - p v + — ) (Rq - a 6 ) . (14) 



Unfortunately if a ^ then there is no simple analytical expression for the critical Blake 
radius R c = R c (a, n,a). In the isothermal case the critical radius can be found explicitly as 
the solution to the cubic equation 
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It is straightforward to show from the discriminant of the cubic equation when a > and G > 
that the equation will have one real solution R c and a pair of ignorable complex conjugate 
solutions. In the case of no hardcore radius, i.e. a = 0, then 
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so that the critical value for the liquid pressure ([13]) is then 



On combining (|15|) and (I16p the Blake critical radius is then given in the familiar form by 
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Note that 
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and similarly 
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so the critical pressure p c and the ambient pressure p^ in general differ by terms O (e) but 
when isothermal O (e 2 ) . Thus when the equilibrium radius is close to the critical radius and the 
ambient pressure is close to the critical pressure, the natural frequency can be expressed as 

2 _ 2 ^ (3* ~ 1) m) 
W ° " p-R% • (20) 



When the perturbation is applied to the Rayleigh-Plesset equation (I10p in this regime both 
the O (1) and O (e) terms are zero if and only if the system is isothermal, i.e. k = 1. Thus, when 
isothermal to O (e 2 ) the governing equation is a Helmholtz oscillator 



N 



x + 2(± + x (1 - x) = A n sin (fi„ (r + r )) , (21) 
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where the overdot represents the derivative with respect to the nondimensional time r and 
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are constant terms of O (1) when e is small for bubbles whose initial radius is of order microns 
driven by frequencies in the megahertz range [10]. When k = 1 the analysis holds to O (e 2 ) 
because thermal dissipation is assumed to be negligible. Here £ is the nondimensional damping, 
O n are harmonics of the nondimensional frequency, A n the nondimensional Fourier components 
of the applied pressure and To is the phase of the incident wave. Note that the series expansion 
is truncated to N terms in order to disregard terms which are greater than O (e 2 ) . 

When forcing and viscosity are rescaled as small parameters, £ h- ► e£ and / i— > ef, which both 
destroy the integrable Hamiltonian structure then the Mel'nikov integral can be calculated. Let 
the e-perturbed system be given by x = f (x) + ef 1 (x, r) with x = (x, y) T = (x, x) T and fi 
an O-periodic function. Explicitly the vector field is given by 

x = y, 



y = x (x - 1) + e ^2 A n sin M - 2 £yj 



The unperturbed system, e = 0, admits a homoclinic orbit 7 emanating from the saddle at (1,0) 
of the form 

l{r) = \ (tanh 2 (I) - l,3tanh sech 2 (I 

The first order Mel'nikov integral at the homoclinic energy level h can simply be calculated 
using Cauchy's residue theory as 
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Due to the summation, it is not possible to formulate an explicit condition for simple zeros, but 
instead perform a calculation numerically. For the Bessel-Fubini and the Fay solutions it is simple 
to show numerically that the Mel'nikov integral has simple zeroes for larger A n than in the case of 
linear wave propagation. A Poincare section can be constructed which is topologically conjugate 
to a Bernoulli shift map on two symbols - in effect the likelihood of a i? max being greater or 
smaller than the previous cycle is as random as the toss of a coin. Thus, for sufficiently small e 
chaotic bubble oscillations will occur in the vicinity of 7. From the simplified normal form, the 
effect of nonlinear wave propagation implies that a cascade to chaos and unbounded bubble 
growth will occur at higher driving pressures P or larger initial radii Rq than for linear wave 



propagation. Note that the classical first order Mel'nikov integral will typically overestimate the 
threshold by not including higher order contributions. 

For the Bessel-Fubini solution as r s —* the linear result is recovered and no summation is 

required, that is as r s — > so P\ — > P, Pj — > for all j = 2, 3, As the distance towards the 

shock decreases, i.e. r s increases from zero, so a higher P or larger Rq is necessary in order to 
have simple zeroes. For the Fay solution higher pressures are required further from the shock. 
For low Gol'dberg numbers, then higher pressures are required than for materials with high 
Gol'dberg numbers. 

As — > oo then A^j, — > 12£/5 and higher order contributions vanish, so that the stable 

and unstable manifolds stay disjoint. However, as — > then A4l — > 12£/5, but now second 
order terms will affect the threshold criteria [18]. Indeed, numerical analysis suggests that higher 
order contributions must be taken into account in this regime. Furthermore, the authors [18] 
state that higher-order Mel'nikov analysis is always necessary when considering nonlinear wave 
propagation due to the interactions of the differing harmonic excitations. 

It is simple to show that in the absence of (nondimensional) viscous damping, i.e. £ = 0, that 
generically the Mel'nikov integral will have simple zeros for all (non-zero) parameter values. The 
effect of viscosity actually reduces the likelihood of violent cavitation, a result also found by 
Szeri [13]. No such analysis has yet been performed in the visco-elastic case but it is currently 
under investigation. 

The bifurcation diagrams presented in figures [3] and [U were computed for the Rayleigh-Plesset 
equation (flQl) with an adaptive, explicit Runge-Kutta method of order (8, 5) due to Dormand 
and Prince [19]. In order to ignore transient behaviour the first fifty cycles were disregarded and 
the next fifty cycles plotted. Figure [3] generalises the inferences from the radius-time profiles 



displayed in figure EJ Subfigure 3(a) illustrates that the maximum amplitude of the bubbles 
under forcing from nonlinear wave propagation is smaller than the maximum amplitude of 
bubbles under forcing from linear wave propagation. By defining £ as the phase of the minimum 
radius R m in = R (imin) per acoustic cycle [12] 

£ — (^min i) fi 

subfigure [3(b)] illustrates that inertial cavitation occurs earlier in each acoustic cycle for nonlinear 
waves than for linear waves. Both bifurcation diagrams in figure [3] show the same regions of 
developing instability as both solution measures determine Poincare sections tangential to the 
motion of the cavity, i.e. when R = 0. 

Figure S] shows a bifurcation diagram beyond the critical threshold and figure [5] then 
shows a selection of associated radius-time profiles illustrating the cascade to chaos through 



a period doubling bifurcation. Subfigure 5(a) shows a stable one-period radius-time profile 



at P = 1.30atm, which undergoes a period doubling bifurcation so that at P = 1.35atm in 
subfigure |5(b)| the profile is two-periodic. For P = 1.375atm as subfigure 5(c) illustrates the 



radius-time profile is four-periodic. By P = 1.40atm the pressure has passed an accumulation 
point of period doubling bifurcations and as subfigure |5(d)| illustrates the bubble oscillations 
are chaotic. Note that the pressure intervals between periodic cycles decreases as the system 
approaches chaotic oscillations, characteristic of the cascade to chaos. From a clinical perspective, 
when the duty cycle maybe in the order of seconds rather than micro seconds, accurate 
predictions for the heat deposition on tumour sites over all but the shortest time scales can 
not be inferred from cavitation activity when the driving pressure is beyond the threshold. 

Bifurcation diagrams of linear and nonlinear wave propagation at equal driving pressures far 
beyond the threshold are significantly different, with differing regions of stability, indicating that 
beyond the threshold predictions based on linear wave propagation will be inaccurate. 



5. Conclusion 

In this paper the effect of nonlinear wave propagation has been investigated numerically and 
analytically. The effect of nonlinear wave propagation redistributes energy from the primary 
harmonic to higher harmonics. From this two significant conclusions can be drawn; firstly the 



maximum bubble radius is reduced. This is clearly illustrated by the bifurcation diagram 3(a) 



Thus the likelihood of shape instability is reduced. Secondly, inertial cavitation occurs earlier 
in each cycle, again as illustrated by the bifurcation diagram |3(b)| The earlier onset of collapse 
allows for the bubble to return to an equilibrium radius and for inertial cavitation to reoccur. 

Amongst the many threshold criteria applied to cavitation, it is worthwhile emphasising, 
although the Mel'nikov criteria does not exactly correspond to the escape boundary whereby 
bubbles increase without bound [20], it initiates the penetration of escaping tongues into the safe 
basin and therefore constitutes the first event in a well known sequence of complicated events 
which leads to unbounded bubble growth. Thus the Mel'nikov criteria provides a lower bound 
threshold value. In the context of therapeutic applications, where safety is paramount, such a 
criteria is of great value. 

In therapeutic applications a cloud of bubbles will exist, comprised of many thousands of 
bubbles of differing equilibrium radii and resonance frequencies. Thus it could be assumed that 
stability criterion for a single bubble is of limited use. However, the stability of a single bubble 
does provide insight into the stability of entire bubble clouds: experimental evidence [21] suggests 
that entire bubble clouds undergo period doubling cascade to chaos, not just the subset of 
bubbles which satisfy the Mel'nikov criterion. It is believed that the interaction between the 
bubbles results in an averaged behaviour. Indeed, numerical and experimental calculations of 
the dimension of the attractor in phase space are remarkably low, between 2 and 2.5, indicating 
that the number of relevant degrees of freedom in the system is also low [22,23]. 

In many applications of therapeutic ultrasound the tissue will have distributed 
inhomogeneities, leading to dissipation through diffraction, but it is conjectured that material 
nonlinearity of is greater significance than material inhomogeneity [24], although this subject is 
currently under investigation. 
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Figure 3. Bifurcation diagram for a bubble determined by the Rayleigh-Plesset equation (|10|) 
with parameters given in figure[2]but initial radius i?o = 0.9/im. The figure shows the normalised 
maximum amplitudes for fifty cycles after fifty cycles in order to disregard transient behaviour. 
The first ten terms of the Bessel-Fubini solution were computed. The radii from linear wave 
propagation are shown in red, those from nonlinear wave propagation shown in blue. 
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Figure 4. Bifurcation diagram showing the period doubling and the onset of quasi-periodic 
oscillations as the magnitude of the forcing pressure of the Bessel-Fubini solution @ is varied 
for the Rayleigh-Plesset equation (jlOp . The parameters are the same as figure [3] except that 
now i?o = 1.4^m. Profiles of the radius time curves at the marked pressures are illustrated in 
figure [5j 
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Figure 5. Period-doubling cascade to chaos as illustrated by the four subfigures marked 



on the bifurcation diagram HJ Subfigure 5(a) shows a stable one-period radius-time profile 
at P = 1.30atm, which undergoes a period doubling bifurcation so that at P = 1.35atm in 
subfigure 5(b) the profile is two-period. For P = 1.375atm as subfigure 5(c) is four period before 
shows a chaotic profile at P = 1.40atm. 
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